Please ignore until reviewed by @ChrisRackauckas.

Follow-up to #236 / Steven's numerical-stability comment.

Summary

• Computes the backward-error scale as exp((logcoeff - delta*log(normTs))/K) instead of (_phi_be_coeff(m, p) / normTs^delta)^(1/K).
• Avoids spurious Inf, NaN, or zero intermediates for large p / tiny normalized norms when the final rooted scale is representable.
• Adds BigFloat-referenced regression coverage for the scale regimes that reproduced the issue (normTs = 1e-8, 1e-2, 0.5, 1.0, 1.9; p = 60, 200, 500).

Local verification

• timeout 3600 /home/crackauc/.juliaup/bin/julia +1.10 --project=. -e 'using Pkg; Pkg.test()'
  • Core/basictests.jl | 673 pass, 1 broken, 674 total
• timeout 3600 /home/crackauc/.juliaup/bin/julia +1.10 -e 'using Pkg; Pkg.activate(; temp=true); Pkg.add("Runic"); using Runic; exit(Runic.main(["--check", "."]))'
  • passed
• timeout 3600 /home/crackauc/.juliaup/bin/julia +1.10 --project=test/qa -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate(); include("test/qa/qa.jl")'
  • Quality Assurance | 15 pass, 1 broken, 16 total

Extra stress checks

Ran a separate local harness against scalar-series / legacy / block-exponential references for:

• tiny diagonal matrix, p = 80: worst relerr 2.78e-17
• tiny dense matrix, p = 60: worst relerr 2.22e-16
• defective Jordan small eigenvalue, p = 8: worst relerr 6.32e-16
• oscillatory complex matrix, p = 6: worst relerr 3.34e-14
• large negative spectrum, p = 10: worst relerr 6.66e-16
• random scaled nonnormal matrix, p = 6: worst relerr 2.68e-15
• workspace high-p call returned cache.info == 0 and finite outputs.